function stats = checkibp(N, alpha, samplerid)
tEKp = alpha*hrmsum(N);
stats.tEKp = tEKp;
% http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm
EKp = 0;
M2Kp = 0; 
Enzr = zeros(N,1);
M2nzr = zeros(N,1);

if samplerid == 0
    % basic sampler
    nsmpl = round(N*tEKp)*10;
    for i = 1:nsmpl
        Z = ibprnd(N, alpha);

        xKp = size(Z,2);
        deltaKp = xKp-EKp;
        EKp = EKp + deltaKp/i;
        M2Kp = M2Kp + deltaKp*(xKp-EKp);

        xnzr = sum(Z,2);
        deltanzr = xnzr-Enzr;
        Enzr = Enzr + deltanzr/i;
        M2nzr = M2nzr + deltanzr.*(xnzr-Enzr);
    end
    VKp = M2Kp/(nsmpl-1);
    Vnzr = M2nzr/(nsmpl-1);

    stats.EKp = EKp;
    stats.VKp = VKp;
    stats.Enzr = Enzr;
    stats.Vnzr = Vnzr;
elseif samplerid == 1
    % conditional sampler
    Z = ibprnd(N, alpha);
    nsmpl = round(N*tEKp)*10;
    for i = 1:nsmpl
        Z = ibpcrnd(Z,alpha,randperm(N));

        xKp = size(Z,2);
        deltaKp = xKp-EKp;
        EKp = EKp + deltaKp/i;
        M2Kp = M2Kp + deltaKp*(xKp-EKp);

        xnzr = sum(Z,2);
        deltanzr = xnzr-Enzr;
        Enzr = Enzr + deltanzr/i;
        M2nzr = M2nzr + deltanzr.*(xnzr-Enzr);
    end
    VKp = M2Kp/(nsmpl-1);
    Vnzr = M2nzr/(nsmpl-1);

    stats.cEKp = EKp;
    stats.cVKp = VKp;
    stats.cEnzr = Enzr;
    stats.cVnzr = Vnzr;
else
    error('Unknown samplerid!');
end

end
